Modeling suggests that virion production cycles within individual cells is key to understanding acute hepatitis B virus infection kinetics

Hepatitis B virus (HBV) infection kinetics in immunodeficient mice reconstituted with humanized livers from inoculation to steady state is highly dynamic despite the absence of an adaptive immune response. To recapitulate the multiphasic viral kinetic patterns, we developed an agent-based model that includes intracellular virion production cycles reflecting the cyclic nature of each individual virus lifecycle. The model fits the data well predicting an increase in production cycles initially starting with a long production cycle of 1 virion per 20 hours that gradually reaches 1 virion per hour after approximately 3–4 days before virion production increases dramatically to reach to a steady state rate of 4 virions per hour per cell. Together, modeling suggests that it is the cyclic nature of the virus lifecycle combined with an initial slow but increasing rate of HBV production from each cell that plays a role in generating the observed multiphasic HBV kinetic patterns in humanized mice.


Introduction
Despite the availability of an effective vaccine for the prevention of hepatitis B virus (HBV), HBV infection continues to impose an enormous burden with an estimated of 270 million chronically infected individuals and about 1 million deaths every year due to complications of HBV, including cirrhosis and liver cancer [1]. Research to elucidate the molecular mechanisms that regulate the HBV lifecycle and infection outcome (i.e., clearance vs. persistence) has been hampered by the lack of model systems that recapitulate HBV infection [2]. Significant attempts have been made to develop small animal models of HBV infection [3]. The most successful small animal HBV infection model approach is based on liver-repopulation in immunodeficient mice repopulated with primary human hepatocytes as these are the natural target of HBV [4][5][6][7][8].
We previously assessed acute HBV infection kinetics from infection initiation to viral steady state in 42 chimeric urokinase-type plasminogen activator transgenic/severe combined immunodeficient (uPA-SCID) mice reconstituted with human hepatocytes [9]. Serum HBV DNA was measured at varying intervals starting as early as 1 min post-inoculation and going out 63 days [9]. Despite varying HBV doses (10 4 −10 8 genome equivalents) and different batches of human hepatocytes, a consistent pattern of distinct phases was observed [9] (Fig 1). To elucidate the processes that result in the complex HBV kinetics observed from initiation to steady state in these mice, we have developed an agent-based modeling (ABM) approach that considers the cyclic nature of the viral lifecycle within individual cells and the resulting distinct waves of viral release and spread. Specifically, the incorporation of viral production cycles within individual infected cells starting with an initially slow production (1 virion every 20 hours) that increases over time to reach steady state of 4 virions every hour recapitulates the multiphasic kinetic patterns. Using ABM to more accurate conceptualize HBV production as cycles rather than a continuous increase thus allows us to reproduce the observed HBV infection dynamics in vivo and provides a new modeling approach for simulating viral dynamics during acute infection.

Agent-based modeling of HBV dynamics
While we showed [9] that the standard experimental approach of measuring HBV DNA levels in serial serum samples from infected mice allows for the calculation of average viral parameters and demonstrates complex multiphasic pattern of viral amplification from inoculation to steady state (Fig 1), it masks the asynchronous infection of individual cells as the virus spreads [10]. We developed an ABM to investigate the dynamics underlying the observed multiphasic HBV kinetic pattern in humanized mice. The ABM accounts for two types of agents: human hepatocytes (cells) and virus in blood (Fig 2A). The cell agents, which are characterized by their infection stage are represented by a square lattice of 3×10 8 cells i.e., the estimated number of human hepatocytes in the humanized mice [11]. The cell agents can be in one of the following three discrete states: uninfected susceptible target (T), infected cell in eclipse phase (I E ) (i.e., not yet releasing progeny virus), or productively infected cell secreting progeny virus (I P ). Due to the absence of an adaptive immune response in these mice, cell death is not considered in the model, hence the cell number was kept constant throughout the simulation period. The virus agent, V, represents the amount of HBV in the blood and is characterized in the ABM as a single global agent which may infect susceptible target cells with rate constant β to become infected cells in eclipse phase (I E ) or to be cleared from blood with rate constant c (Fig 2A). The ABM execution is an iterative process where each iteration represents a "tick" or a discrete time step, where 1 step = 1 hour.
Because HBV is a noncytolytic virus [12] intracellular viral production increases over time until reaching a resource restriction plateau. To model this on a per cell basis, we quantify both the amount of viral production by infected cells at a given time and the production cycle, i.e., the time interval over which the cell produces virus (Fig 2B). To capture such dynamics, we formulated the following equations.
The amount of virion produced, P, at production cycle τ is determined as: where, P(τ) is number of virions produced by infected cells a τ, P st is steady state virus production, α is number of cycles to reach to 50% of P st , γ is steepness of the production curve, and τ is the production cycle. Viral production, P st , is estimated at steady state when all target cells are infected, therefore P st = cV st /I p where V st represents viral load at steady state, and c represents the virus clearance rate constant from blood. And the interval between cycles is determine by an exponential decay function with: Where l τ is interval between production cycle, τ is the production cycle, δ is scaling factor indicating the initial production cycle length, and ω is decay constant. In the model, both P(τ) and l τ are roundup to the nearest integer. The combination of Eqs 1 and 2 allows for each productively infected cell to slowly produce virions once I E becomes I P with increasing numbers of secreted virions within shorter time intervals until the cell reaches a steady state production.

ABM reproduces multiphasic HBV kinetics after a high dose viral inoculation
The model reproduces well the complex HBV DNA serum patterns observed in 4 mice (Table A in S1 Text) inoculated with 10 8 HBV genome equivalents and followed from inoculation until day 51 post-inoculation (p.i.) (Fig 3 and   , not yet releasing virions), I P = productively HBV-infected cells (i.e., actively releasing virions). The free virus in blood, V, is composed of infectious and non-infections virions. The parameter ρ represents the fraction of virions that are infectious, β represents the infection rate constant, O represents eclipse phase duration, c, represents viral clearance from blood and P(τ) (Eq 1) represents virion secretion from I P . (B) Schematic diagram of viral production cycle for an individual infected human hepatocyte. P(τ) is the number of virions produced by an infected cell, and l(τ)) is the time interval between production cycle (h). The virions were initially released by I P starting with a long production cycle of 1 virion per cell (Stage 1:~0-2.5 days) that gradually reaches a production of 2 virions per cell with a shorten production cycle (Stage 2: 2.5-3 days) and then proceeds to 3 virions per cell (Stage 3:~3-4 days) before virion production increases to reach to a steady state production rate of 4 virions per hour per cell (Stage 4:~4 days onwards). https://doi.org/10.1371/journal.pcbi.1011309.g002

PLOS COMPUTATIONAL BIOLOGY
pictures of the cell populations at distinct time points (Fig 4B) during the simulation of infection for mouse 1 (M1) reveals that while serum HBV DNA is multiphasic (Fig 3), the total cellstate kinetics are simple, e.g., uninfected productively infected cells follow simple sigmoidallike patterns (Fig 4A, green and red curves).

Virus parameter estimates of individual cells that allow for simulation of the observed complex serum viral patterns
The model provides insights into early virus-host dynamics of individual cells from infection initiation to steady state that unmask the nature of the observed complex serum viral kinetic patterns. To fit the data, the model predicts a variable eclipse phase ranging from 5-50 hours in the 4 mice inoculated with 10 8 HBV genome equivalents (Table 1). Using M1 as a representative example, the simulated virion production with two eclipse phase durations of 9 hr (Fig 5A, shaded box) and 48 hr ( Fig 5C, shaded box) illustrates the delay in virion production that results from the longer eclipse phase without affecting subsequent virus production cycles. Post-eclipse phase, the model predicts viral release from productively infected cells starts slowly with a long production cycle of 1 virion per 20 hours that gradually reaches 1 virion per hour (Fig 5B and 5D) after~3-4 days before virion production increases to reach to a steady state production rate of 4 virions per hour per cell (Fig 5A and 5C). A similar picture was found for mice M2, M3 and M4 (Fig V in S1 Text).  Table 1. https://doi.org/10.1371/journal.pcbi.1011309.g003

Dissecting the nature of each serum HBV DNA kinetic phase
Focusing on the details of the model simulation for representative M1, the estimated serum HBV DNA is shown on a time scale that allows visualization of each serum HBV DNA kinetic phase (Fig 6A and 6B). The number of non-productive (I E ) and productive (I P ) infected cells along with the average virion production per productive cells are plotted on the same scale for direct comparison (Fig 6C and 6D).
During the first 6 hours p.i. (Fig 6A, Phase 1), as the model recapitulates the rapid serum HBV DNA (V) clearance from the blood (t 1/2 = 1 hr), it predicts that about 1×10 5 cells were infected, i.e., the first wave of infection. This consists of an initial peak of eclipse phase cells (I E) ) ( Fig  The observed intermediary serum HBV DNA steady state (or plateau) is recapitulated by the model (Fig 6A, Phase 4) as the majority of I P cells are at steady state levels of viral production and the new increasing numbers of I P cells are in the early low level virus production phase of infection (Fig 6C, increasing solid red line). Likewise, the increasing number of I E cells do not contribute to virion production ( Fig 6C, increasing solid gold line) resulting in an overall decrease in the per cell virion production rate (Fig 6C, decreasing green line).
As the infection becomes less synchronized due to the stochasticity that exists in the timing of individual infection events and target cells become limiting, the distinct cycles of infection become less discernable and all subsequent amplification appears as a single viral exponential Table 1. Best estimated ABM parameters for calibrated mice and their estimated space (J<0.7). Parameters are defined in Eqs 1 and 2; GE, genome equivalent; Jscore, represents the objective function (Eq 3) score of the genetic algorithm (GA) fits where the min J-score is the best fit curves shown in Figs 3 and 7. Parameter estimates shown in () represent the min and max of all GA fits with J-score < 0.7 shown in Figs A-G in S1 Text with their statistics in Table C Table 1). The spatial structure on the 2D lattice was not considered in the ABM to affect viral infection spread. https://doi.org/10.1371/journal.pcbi.1011309.g004

PLOS COMPUTATIONAL BIOLOGY
expansion from day~8 until~30 days p.i. (Fig 6B, Phase 5) during which the final target cells are infected and become productive (Fig 6D, solid gold and red lines, respectively) and subsequently progress towards maximal average viral production (Fig 6D, green line). Once all the target cells are productively infected (Fig 6B, solid red line) and achieve maximal average viral production (Fig 6D, green line), serum HBV levels attained steady state (Fig 6B, Phase 6).

ABM reproduces multiphasic HBV kinetics after low dose infection
We previously showed that lower inoculation of 10 7 , 10 6 and 10 4 HBV genome equivalents in humanized uPA-SCID mice led to a delayed, but similar complex HBV kinetic pattern [9]. Therefore, we validated the model against this kinetic data (Table A in S1 Text) by changing only the inoculation dose in the model simulation. Importantly, the model reproduced well the kinetic patterns within the same parameter space estimated for mice that were inoculated with 10 8 HBV genome equivalents (Fig 7, and Figs E-G in S1 Text and Table 1).

Short virion production cycle and eclipse phase lengths diminish multiphasic kinetic pattern
Having validated the ability of the model to accurately recapitulate the kinetics of HBV infection, we proceed to investigate how the two key features of the model, namely cellular eclipse phase and/or production cycle, might affect acute viral infection kinetics in general. To test the effect of shortening and/or lengthening the cellular eclipse phase and/or production cycle, in Lengthening the cellular eclipse phase resulted in an initial period of no viral production followed by a delayed but multiphasic viral amplification (Phases 2-6) (Fig 8A, red dotted  line). In contrast, a short eclipse phase allowed for visualization of more extreme cycling during Phase 2 of the infection without affecting subsequent viral phases (Fig 8A, green dashed  line). Increasing the virion production essentially eliminated the Phase 2 initial plateau resulting in earlier amplification (Fig 8B, green dashed line). Slower production in contrast lengthened the initial lower plateau Phase 2 and delayed subsequent viral phases (Fig 8B, red dotted  line). Interestingly, combining a shorter cellular eclipse with a faster production cycle produces an almost single amplification phase analogous to that observed for many acute viral infections (Fig 8C, green dashed line) perhaps suggesting why such multiphasic viral infection kinetic patterns have not been universally observed when monitoring acute infection kinetics of other viruses. In contrast, combining a shorter cellular eclipse with a slower production cycle produces again the short eclipse extreme cycling this time during a longer Phase 2 followed by delay of all subsequent viral phases (Fig 8C, red dotted line). Combining a longer cellular eclipse with a faster production cycle resulted in an initial period of no viral production followed by immediate amplification (i.e., no Phase 2 initial plateau) with little effect on subsequent viral phases (Phases 3-6) (Fig 8C, green dashed line). While the longer cellular eclipse coupled with a slower production cycle resulted in an initial period of no viral production followed by delayed multiphasic viral amplification (Phases 2-6) (Fig 8D, red dotted line).

Discussion
Standard population-based measurements of viral infection time courses typically show an initial eclipse phase followed by exponential increase ending either in cell lysis or steady state viral levels [13][14][15][16][17][18][19][20][21]. However, we recently reported a surprising complex HBV infection kinetics in chimeric uPA-SCID mice with humanized livers [9]. To investigate the dynamics underlying the unexpectedly complex HBV infection kinetics from inoculation to steady state in humanized uPA-SCID mice, we developed an ABM approach to explain the serum population kinetics from the perspective of individual cell infections. This allowed for the simulation of HBV production from individual cells as cycles rather than a continuous increase which The majority of mathematical models of acute viral infection are based on the assumption of well-mixed virus and cell populations [13][14][15][16][17][18][19][20]. For a noncytolytic virus such as HBV, in the absence of the development of immune response in which infected cells are targeted for loss/ death at a faster rate compared to uninfected cells, these models predict a roughly monophasic viral increase that reaches a high viral load steady state (or peak) once all target cells are infected. Some models also accounted for an eclipse phase, e.g., [22,23], in which newly infected cells remain in a latent phase before becoming virion producing, but even with such additions a viral increase is predicted until all target cells are infected. For example, in the chimpanzee model of acute hepatitis C virus (HCV) infection, viral levels increased in a biphasic manner with a transient viral decline in between (1 week p.i.) concomitantly with the induction of type I interferon [16]. However, HBV is often referred to as a stealth virus because it does not induce significant innate immune signaling [12], making immune signaling a less likely explanation for the shift from rapid expansion (Fig 1, Phase 3) to a slower interim plateau (Fig 1, Phase 4). While intrahepatic interferon induction cannot be ruled out without further experiments, we show here that this shift as well as all the phases observed can be accurately recapitulated by describing viral dynamics at the individual cell level combining the viral eclipse phase, the increasing rate of virus secretion based on increasing production cycles, and the resulting waves of new infection (Fig 6).
Viral kinetics early post-infection often exhibits a viral decline phase followed by lower plateau or undetectable viral level (a.k.a. eclipse phase) before exponential amplification. This is assumed to reflect the time period in which clearance of input virus in the blood is balanced by de novo viral production. We showed that mice inoculated with high (10 8 ) virus inoculation HBV DNA had a low viral plateau which lasted approximately 1-3 days p.i (Fig 1, Phase 2). Fitting the ABM with measured serum HBV DNA we estimated a range of 5 and 50 hours in which newly infected cells (I E ) cell remain in a latent phase before becoming virion producing infected cells (I P ). While this initially seemed like an unexpectedly broad parameter range, single cell analysis of viral infections has revealed analogously large ranges of variability in the progression of viral replication [24,25].
We previously reported [9] that HBV DNA inoculum size had no effect on initial HBV clearance (Fig 1, Phase 1), the viral doubling time during the HBV expansion phases (Fig 1,  Phases 3 and 5), the length of the interim plateau (Fig 1, Phase 4), or viral steady-state levels (Fig 1, Phase 6), but rather resulted in a lower viral plateau (Phase 2) and a delay in detection of initial virus expansion (Fig 1, Phase 3), which subsequently delayed all other kinetic phases. Importantly, the model confirms this as simulating the infection kinetics in mice after HBV inoculation of 10 7 , 10 6 or 10 4 genome equivalents, requires no significant change in any of the estimated ABM parameters compared to the mice inoculated with 10 8 genome equivalents except for the inoculation dose itself (Table 1), further supporting our ABM modeling approach.
Finding that the unusual multiphasic kinetics observed can be reproduced by a model that is based on the inherent cyclic nature of a viral lifecycle, raised the question why such complex kinetics were observed for HBV but not more broadly for other viral infections. While it is possible that such a multiphasic pattern could simply be missed in the absence of frequent sampling, this complex pattern was also eliminated by increasing viral production and reducing the eclipse phase in the model simulations ( Fig 8C). Notably, these parameter changes are consistent with the faster lifecycles associate with many acute viruses and may explain why multiphasic viral kinetics has not been routinely observed for other viruses.
The current ABM does not account for intracellular HBV-host dynamics that explain the estimated timing of the viral production cycles but provides a prediction that can now be investigated. Our previous kinetic study [9], indicated that intrahepatic total HBV DNA, cccDNA, and RNA correlate with serum HBV DNA during infection in these mice. Thus, a future detailed intracellular kinetic analysis may shed light whether intrahepatic HBV RNA and/or DNA levels exhibit the same multiphasic amplification pattern and how cccDNA recycling impacts the level of HBV secretion early in infection.
Meanwhile, we show in the current study that the incorporation of viral production cycles into an ABM recapitulates the multiphasic serum HBV kinetic patterns observed in uPA-SCID mice reconstituted with humanized livers from inoculation to steady state [9]. Importantly, such complex HBV infection kinetics can be seen also in immunocompetent chimpanzees [14,26] indicating that this complex picture is not unique to humanized uPA-SCID mice. Thus, using agent-based modeling to more accurate conceptualize virus production as cycles rather than a continuous increase allows us to reproduce the observed HBV infection dynamics in vivo and provides a new computational approach for simulating viral dynamics during acute infection.

Ethics statement
All animal protocols from which the data in this manuscript were derived were performed in accordance with the Guide for the Care and Use of Laboratory Animals and approved by the Animal Welfare Committee of Phoenix Bio Co., Ltd (approval number 1031).

Mice and HBV
The in vivo data investigated in silico herein has been previously published and described in detail in [9]. Briefly, humanized liver chimeric mice were produced by splenic injection of cryopreserved human hepatocytes into uPA/SCID mice as previously described [8,9]. Human hepatocyte repopulation rates were estimated by blood human albumin levels and chimeric mice showing repopulation rates of greater than 70% were injected with serum containing 10 4 , 10 6 , 10 7 and 10 8 copies HBV DNA (genotype C, Accession No. AB246345), originally provided by Dr. Sugiyama [27] via the tail vein. Preparation of the inoculum and methods utilized for DNA extraction and quantification of HBV DNA were performed as previously reported [9,28].

Parameter estimations
We previously showed that during the first 6 h p.i., and 6-24 h p.i. HBV was cleared from blood with a t 1/2 of~1 h and~3 h, respectively, independent of inoculum size which ranged from 10 4 to 10 8 cp/ml [9]. Because the above two initial clearance phases have been combined into one, now jointly termed Phase 1 (Fig 1), we assumed a fixed clearance rate of virus from blood as c = 0.5 h -1 . The fraction of virus in the blood (V) that was infectious was arbitrary set as ρ = 0.5. Based on the experimental data, HBV serum levels at steady state were set at, V st = 9.3±0.3. Viral production, P st , is estimated at steady state where all target cells are infected (I p = 3×10 8 cells), which is equivalent to P st ¼ 0:5X10 9:3�0:3 3X10 8 � between 2 and 7 virion/cell. The remaining parameters were estimated by calibrating the ABM with the experimental data (Table 1), using parameter constraints as shown in Table B in S1 Text that were obtained by preliminary ABM simulations and fitting with the experimental data using Anylogic.

Model calibration
Model parameter fitting was done using a Genetic Algorithm (GA) [29,30] with the EMEWS framework [31] on the Midway2 high-performance computing (HPC) cluster at the University of Chicago. Midway2 has 400 nodes, each with 28 cores and 64GB of memory. Some additional development was done on the Bebop HPC cluster, managed by the Laboratory Computing Resource Center at Argonne National Laboratory. Bebop has 1024 nodes comprised of 672 Intel Broadwell processors with 36 cores per node and 128 GB of RAM and 372 Intel Knights Landing processors with 64 cores per node and 96 GB of RAM. The GA was implemented using the DEAP [32] evolutionary computation Python framework (specifically [33]: Chapter 7) and integrated into an EMEWS HPC workflow using EMEWS queues for Python (EQ/Py) [31]. The input parameters for the calibration were all in linear space. The objective function used a log transform of the model outputs: J ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 n X ½log 10 Vðt i Þ À log 10V ðt i Þ� where log 10 V(t i ) are the log-transformed predicted values of the viral load, log 10V ðt i Þ are the log-transformed experimental values, σ(t i ) are the standard errors of the experimental data, which are assumed to be 0.5 for each data point t i , and n are the number of sample data points. The use of HPC resources enable the concurrent evaluation of large numbers of design points (n = 102), reducing the time to solution. During each iteration of the GA, the best points from the currently evaluated population are selected using a tournament selection method to create a new population. Each of these points is then mated with another according to a crossover probability and, finally, each of the resulting points is mutated according to a mutation probability. At each GA algorithm iteration, the new population is evaluated in parallel, and the evaluation results are gathered. The GA population size was set to 102, the mutation probability to 0.2, the crossover probability to 0.5, and the number of iterations to 25. On Midway2 the runtime for a typical run was 7.3 hours using full concurrency on 28 nodes (with 28 cores per node), or about 5700 core hours.